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We study the effects of a phase difference on Yu-Shiba-Rusinov (YSR) states in a spinful Coulomb- 
blockaded quantum dot contacted by a superconducting loop. In the limit where charging energy 
is larger than the superconducting gap, we determine the subgap excitation spectrum, the cor¬ 
responding supercurrent, and the differential conductance as measured by a normal-metal tunnel 
probe. In absence of a phase difference only one linear combination of the superconductor lead 
electrons couples to the spin, which gives a single YSR state. With finite phase difference, however, 
it is effectively a two-channel scattering problem and therefore an additional state emerges from the 
gap edge. The energy of the phase-dependent YSR states depend on the gate voltage and one state 
can cross zero energy twice inside the valley with odd occupancy. These crossings are shifted by the 
phase difference towards the charge degeneracy points, corresponding to larger exchange couplings. 
Moreover, the zero-energy crossings give rise to resonant peaks in the differential conductance with 
magnitude 4e^//i. Finally, we demonstrate that the quantum fluctuations of the dot spin do not 
alter qualitatively any of the results. 

PACS numbers: 72.10.Fk, 74.45.+C, 73.63.Kv, 74.50.+r 


I. INTRODUCTION 

Yu-Shiba-Rusino\ffH3l states are bound states in an s- 
wave superconductor induced inside the energy gap by 
local magnetic moments. Individual localized Yu-Shiba- 
Rusinov (YSR) states have been observed both by scan¬ 
ning tunneling spectroscopy of magnetic atoms like Mn or 
Gr adsorbed on superconducting Pb or Nb substrates,®^!! 
and as subgap states in Goulomb blockaded ^antum 
dots (QD) coupled to superconducting (S) leads.®^The 
quantum dot realization is based on the spin-1/2 of odd- 
occupation charge states, and is therefore free of most 
of the material dependent complications for adatoms on 
a surface, like mixed valence, higher spin, and magnetic 
anisotropy. Furthermore, the quantum dot system al¬ 
lows for electrical tunability of the particle-hole asymme¬ 
try and, to some extent, the exchange coupling between 
the spin on the quantum dot and the quasiparticles in 
the superconductor, which makes it an ideal system for 
studying the properties of individual YSR states. 

Recently, sharp subgap states have been observed with 
a weakly coupled normal (N) metal lead probing an S- 
QD-S junction as illustrated in Fig. where the quan¬ 
tum dot was formed in a Goulomb-blockaded segment of 
carbon nanotub^^^or InAs nanowir^I^ spanning a super¬ 
conducting A1 loop. The gate control allows for determi¬ 
nation of even (spinless) or odd (spinful) charge states 
of the quantum dot by even-odd effects (such as absence 
or presence of Kondo resonance^^^ of the stability dia¬ 
gram. By tuning the magnetic flux piercing the A1 loop, 
one may control the phase-difference across the quan¬ 
tum dot, and this device therefore provides additional 
information about the phase dependence of the subgap 



Fig. 1. (Color online) Sketch of the device studied here (fol¬ 
lowing the experiments in Refs. [12] and El, comprised of a 
Coulomb-blockaded QD (green), tunnel coupled to supercon¬ 
ducting leads (blue) with an applied phase difference, 0, con¬ 
trolled by a magnetic flux. A backgate voltage is adjusted to 
provide a single spin-1/2 on the QD, and the resulting phase- 
dependent YSR states are probed by a normal-metal tunnel 
probe (yellow). 


states. It is this phase dependence of the spin-induced 
YSR states, which is addressed in the present paper. 

Earlier measurements of super cur rent thr oug h phase- 
biased Coulomb-blockaded quantum dots,^^^^^ have al¬ 
ready demonstrated that odd-occupied spin-degenerate 
dots may lead to negative (7r-phase) super current. Also, 
a supercurrent sign reversal, i.e., a tt — 0 transition, has 
been shown to take place when adjusting the gate voltage 
to move away from odd occupancy, thereby increasing the 
ratio of the Kondo temperature Tk to the superconduct¬ 
ing gap A, in accordance with a number of theoretical 
predict ions Whereas the supercurrent only carries 
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information about the ground state, the additional spec¬ 
troscopic information from experiments like Refs. [12] and 
[T3l now offers an opportunity to learn more about the 
subgap excitations caused by a local magnetic moment. 

The nature of subgap states depends on the ratio of 
the charging energy Ec and the superconducting gap A. 
For Ec > A and Ec ^ F, where F denotes the elastic 
broadening of the dot states due to their coupling to the 
leads, the natural starting point for the spinful dot is the 
Kondo modeP^ with exchange coupling J and with the 
lead conduction electrons described by BCS Hamiltoni¬ 
ans. As we show below, in absence of a phase difference 
across the junction, this reduces to a single-channel YSR 
problem. For weak exchange coupling, Tk ^ A, where 
the ground state is a doublet, the YSR state is an ex¬ 
cited singlet state consisting of a single quasiparticle in 
the lead bound to the dot spin. As the exchange coupling 
is increased the excited state crosses zero energy and 
the ground state changes abruptly to spin singlefP^^®] at 
Tk ^ A, where Tk is the Kondo temperature. For even 
larger values Tk ^ A, the ground state can be described 
as the well-known Kondo singlet. 

In the opposite limit, Ec < A, the natural starting 
point for understanding the subgap states is a model 
where the superconducting electrons are integrated out, 
giving rise to a local pairing on the dot with strength 
F, which was studied by Meng et alW^ It gives a mixing 
of states with occupation N ± 1 when the average oc¬ 
cupation is N. With odd average occupation there will 
be two subgap states split by F, as found already in the 
single orbital Anderson model.[^ Even though the two 
cases Ec < A and Ec > A are naturally described in 
different languages, the physical situations are similar. 
In both cases, the excited states correspond to an extra 
quasiparticle bound by the local spin. For Ec > A the 
bound particle resides mainly in the superconductor in 
the form of a YSR state, while for Ec < A it resides 
mainly on the dot because of the hybridization of 0 and 
2 electrons. The two situation are illustrated in Fig. 
For the experiments reported in Refs. [13] and nn the rel¬ 
evant limit is A < Ec^ which is also the limit considered 
in this paper. 

For a multi-level quantum dot, the finite level spac¬ 
ing, 5E^ will also affect the simple evolution of ground 
states when 5E ^ Tk, A. In this paper we thus restrict 
our attention to small dots, for which the level spacing 
is the largest energy scale, and F, A ^ Ec ^ SE, com¬ 
plementary to the A ^ Ec regime studied by Meng et 
al^ Furthermore, we restrict our attention to gate volt¬ 
ages adjusted to accommodate odd occupation and hence 
a net spin-1/2 on the quantum dot, which will be de¬ 
scribed within an effective cotunneling (Kondo) model. 
For the main part of the paper, however, the dot spin is 
treated as classical (non-fluctuating), being polarized in 
a fixed direction. We will show that this approximation 
does not modify the physics substantially (at least for 
weak dot-lead coupling), while considerably simplifying 
the calculations. 



Fig. 2. (Color online) Schematics of the two situations, 
Ec > A and Ec < A. In both cases there is a subgap 
state in the odd diamond. For large A it appears because 
of hybridization of the iV = 0 and N — 2 charge states by 
the superconductor, while for large Ec it appears because of 
hybridization between the dot electron and quasiparticles in 
the leads, forming a bound state, which is the YSR state. 
In both cases the average occupation on the dot is equal to 
one and the subgap states have similar dependences on gate 
voltage. In this paper, we focus on the situation (a) with 
Ec > A. 


The paper is organized as follows. In Section [H] an ef¬ 
fective Kondo model for an odd-occupied quantum dot 
coupled to superconducting leads is introduced. In Sec¬ 
tion |TTI| we derive the subgap states within the polarized- 
spin approximation. The corresponding nonlinear tun¬ 
neling conductance to an auxiliary normal lead is dis¬ 
cussed in Section EISL with a few technical details 
relegated to Appendix |A[ and finally the supercurrent 
through the phase-biased S-QD-S junction is addressed 
in Section |III C[ In Section [TVl we briefly compare the 
results from the polarized-spin approximation with the 
perturbative (in dot-lead coupling) results including the 
full quantum dynamics of the dot spin, which we derive 
in Appendix [B] 


II. THE MODEL 

We consider a quantum dot tunnel-coupled to two su¬ 
perconducting leads, capacitively coupled to a gate elec¬ 
trode, and subjected to an applied magnetic field. The 
coupling to the normal metal electrode is also included 
for the purpose of doing transport spectroscopy. 

The quantum dot connected to these three leads is 
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modelled by a single orbital Anderson-type model: 


where the Zeeman term for the quantum dot spin reads 


= + ( 1 ) 

The superconducting leads are described by the effec¬ 
tive Bardeen-Cooper-Schrieffer (BCS) mean-field Hamil¬ 

tonian 

ikanaka^aka (^) 

aka 

~ ^ + ^a^a-kl^akt) ’ 

ak 


(8) 

with S'* denoting the spin operator on the dot.The trans¬ 
formation is valid for A/[/ <C 1 and \gce— gd\B/U <C 1.^ 
For the polarized-spin approximation considered in the 
next section, the Zeeman term for the dot spin has no 
influence, but it will become important in Appendix [B] 
where the magnetic field dependence of the quantum cor¬ 
rections are discussed. 

The exchange cotunneling term reads 


where 4^ = Cfe + (^ 9 ceBl 2 and 4 = s^k - 9a are the 
conduction electron dispersions. Here a = L/R labels 
the two superconducting leads and k and a are lead or¬ 
bital and spin quantum numbers, respectively. The leads 
are assumed to have bandwidth 2D with ^ [—D^D]. 
The Lande ^-factor in the superconductors is denoted by 
the magnetic field by H, and the complex supercon¬ 
ducting order parameters by Aq, = Ae'^^^. Here we for 
simplicity take the two order parameters to be of equal 
magnitude. The normal-lead Hamiltonian is 

= ^(^AT/ccr — /^A^)c]v/ccr%/ccr’ (^) 

ka 

with a bias voltage applied to this lead /in — ha = 
Assuming the level spacing in the quantum dot to be 
large, we model it by a single orbital: 


Hj= Y. Ja'aS^^lk'a'rk.C^ka, (9) 

i=x,y,z 
a' k' a' aka 


where r* denotes the Pauli matrices, and the potential 
scattering term is 

Hw = Y ^a'a^i'k'a^^aka- (10) 

a' k' ,ak,a 


Here the exchange, and potential scattering amplitudes 
are given by 


Jaa' 


tatn 


aa' — 


2X tcyta 


( 11 ) 


l_a;2 u ’ ““ 1-^2 U ’ 

where x parameterizes the dimensionless gate voltage as 


Hd= Y Sddid^ + Un^ni, Sda = £d + (r9dB/2, ( 4 ) 


X = 1 


U ‘ 


( 12 ) 


where dj. creates an electron in the orbital with spin <j, U 
is the charging energy on the quantum dot, Cd is the level 
position, which is controlled by the gate voltage, and g^ 
is the ^-factor in the dot (which can be different from 
dee)- 

The coupling between the leads and the dot is de¬ 
scribed by the tunneling Hamiltonian 

HT=Y{toclkada+t*adic.ka), ( 5 ) 

aka 


Note that the Anderson model always gives antiferromag¬ 
netic exchange J > \W\ >0 inside the odd occupied di¬ 
amond X G [—1,1] and that W breaks particle-hole sym¬ 
metry and therefore vanishes at the particle-hole sym¬ 
metric point, X = 0, defining the middle of the Coulomb 
diamond. We also define the following dimensionless cou¬ 
pling constants: 

Qaa' = TT^pJaa'S, g = gLL + gRR, 

Waa' = TriypWaa' , ^ = ^LL + ^RR, 


where t^ denote the lead-dot tunneling amplitudes and 
the lead index a is extended to run through L, i?, and 

N. 

We focus on the odd-occupied spinful cotunneling 
regime well inside the corresponding Coulomb diamond, 
where 

2TrPF\ta\‘^ =ra -Sd,U + Sd, (6) 

with Up denoting the density of states at the Fermi level. 

A standard Schrieffer-Wolff transfer mat iorP^^ then 
leads to the following effective cotunneling (Kondo) 
model for the spin-^ coupled to the normal and two su¬ 
perconducting leads: 

(7) 


to be used extensively below. 


III. POLARIZED-SPIN APPROXIMATION 

We start by considering the case where the spin op¬ 
erator in Eq. is treated as a classical variable with 
a fixed direction: S ^ Sz. In this approximation, the 
problem is similar to the original problem considered by 
Yu, Shiba, and Rusinov,IIH^ but now with two supercon¬ 
ductors having different phases. In Appendixwe show 
that this approximation is justified when g^B g‘^ A. At 
zero field the excitation energies calculated within this 
approximation correspond to the correct result at weak 
coupling {g 1) only after rescaling ^ by a factor of 


H = i^N + 77 lr + + Hj + Hw, 
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3. For now, we use the polarized-spin approximation in 
order to discuss the spectroscopy. 

We start by diagonalizing the exchange, and potential 
scattering terms in L/i^-space, omitting the coupling to 
the normal lead. This diagonalization is possible because 
they share the same matrix structure in L/i^-lead space 

JoiOi' and 1 Tq/q// bFOo/Q//, 


where 


f cos^ 0 sin 0 cos 0 A 
y sin 0 cos 0 sin^ 0 j ’ 


(15) 


with coupling asymmetry parameterized by an angle 0 
defined by 


(cos 6 »,sin 6 >) = {tL,tR)/t, t = \Jt\ + t\. (16) 

Notice that ^ = 7 r /4 corresponds to symmetric coupling 
th = tR. By means of a gauge transformation, the phases 
of the individual pairing potentials in the contacts, Aq,, 
can be combined to a phase difference, 0 = 
appearing only in the scattering terms via the matrix 0 : 



This matrix has eigenvalues 0 and 1 with corresponding 
eigenvectors: 


Vo = 


— COS 0 


Vl = 


cos 0 

g-i0/2sin6l ) ’ 


(18) 


and a unitary transformation that diagonalizes J^a' and 
Waa' is therefore achieved with 


U = 


cosO 

—cosO 


(19) 


leading to the following transformed cotunneling terms: 


Hj + Hw=Y1 


( 20 ) 


kk' (T 


where the new operators co/e and cik are defined as 


Ik — 


^0/c 


= u 


^Lk 

^Rk 


( 21 ) 


From expression (20) it is evident that the channel cor¬ 


responding to eigenvalue 0 does not couple to the cotun¬ 
neling terms. 

However, with finite phase the two channels are not in¬ 
dependent and therefore the problem remains effectively 
a two-channel problem. In order to see this, we write 
lead Hamiltonians in Nambu space 



ik ■^(TgceBl2 

A 


A 

^k T ^QceB/2 


Cka, 

( 22 ) 


with Nambu 4-spinors defined as 

B'ka {^Lkai^Rkai — . (23) 

Notice that since Hq is diagonal in lead space, a unit 
matrix in lead space is implied on each of the four ma¬ 
trix elements in ( 22 ). After the rotation ( 21 ), which in 


Nambu space reads 




the lead Hamiltonian becomes 


TT ( Cfe + crfi'ce-B/2 AP 

"0 APt -a+(T5ceS/2 

ka ^ 


where P = t/’f't/*, which evaluates to 


Cka, 

(25) 


/ cos^ 6 + sin^ 6 —i sin(26>) sin(0/2) \ . . 

y —i sin(26>) sin(0/2) cos^ ^ + 6“*"^ sin^ / 


From the expression for P, we now see that for 0 7 ^ 0 
the two channels are not independent, and the problem 
therefore remains a genuine two-channel scattering prob¬ 
lem. 


A. Subgap spectrum from T-matrix 


We are interested in understanding the structure of 
the subgap states and therefore we study the T-matrix 
for spin-cr given by 

t7(w) = V; [1 - G^ico - ag,,B/2)V,] , (27) 

with the diagonalized exchange, and potential scattering 
terms extended to Nambu space as 


/crJP + VF 0 0 0 \ 

I 0 0 0 0 1 

“ 0 0 aJS - VF 0 ’ 

V 0 0 0 0/ 

whereby 

Hj + Hw = \Y. 

kk' a 


(28) 


(29) 


The local Green’s function corresponding to (25) is found 
as 




f ik 

AP 


<1 

1 

-ik 

)J 


( 

UJ 

0 

TTUr 


0 

_ 5k A 

UJ 


n -1 






\ icA aA 


aA 

—icA 

UJ 

0 


( 30 ) 


—icA \ 
a*A 
0 

CJ / 
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where the /^-integration was performed assuming a con¬ 
stant density of states vp and assuming that jcjl < A ^ 
D. The lead asymmetry and the phase-difference are en- 


coded in the dimensionless coefficients 


a = cos^ 0 + sin^ 6>, 

(31) 

c = sin(26>) sin(0/2). 

(32) 

The condition for poles in is 


det [1 - Gq„{E - (jgceB/2)Va\ = 0, 

(33) 

and after some algebra one finds two roots for 

each spin. 


corresponding to YSR states at energies 


E±,a = cr5ce-B/2 


crc±A 


y(l + 'U)2 +4g2 L 
+ 2y ± 25^52 


(l + u)(n-xw) 
■11/2 

1 

(34) 


with the following shorthand notation 

X = 1 - sin^(26») sin^(</>/2), 

u = w^ -y, 
c_ =sgn(l + X'w), 


(35) 


c+ = 1. 


Here we have assumed ^ > 0; the corresponding solutions 
for ^ < 0 simply have opposite spins, and hence are given 
by 

In the rest of this section we examine the dependence 
of the subgap states on the coupling strengths and the 
phase difference. For zero phase difference, 0 = 0, and 
antiferromagnetic coupling, ^ > 0, Eq. (34) gives the 
following subgap energies: 


Err = -crA 


1 + - g‘^ 




^ > 0. (36) 


This is the result obtained by Yu, Shiba and Rusinov.l^H^ 
In the case when there is no exchange coupling, i.e., g = 
0, Eq. (34) yields the usual expression for the Andreev 


bound state®® 


E 


0 

±,cr 


TcrAWl - 


0 

2 ’ 


w‘^ sin^(26>) 

1 + 


(37) 


where r denotes the normal state transmission of the 
junction. 

The energies of the four bound states are plotted as 
dashed curves against gate voltage and phase difference 
for different cases in Eig. The color scale refers to the 
differential conductance for the system with an added 
normal lead (see next section). In panels (a) and (b) we 
plot the bound-state energies for zero phase difference 
and weak, and strong coupling to the superconductors, 
respectively. The difference between (a) and (b) can be 



Fig. 3. (Color online) The figure shows the YSR state energies 
V = E±^a given in Eq. ( [34| as dashed lines. The YSR en¬ 
ergies are plotted as a function of dimensionless gate voltage 
X = 1 -h 2edlU in (a)-(d), and as function of phase differ¬ 
ence in (e)-(h). x = 0 corresponds to the middle of the odd 
diamond and x = ±1 are the charge degeneracy points, see 
Eqs. 0 and ( |12| ). The lines are overlayed on co ntou r plots 
of the differential conductance (see Fig. [^and Eq. ( [^ ) using 
gNN = 0.1 and the parametrization of the couplings given in 
Eq. ( |A4[ ). The arrows indicate the spin of the excited states. 
The exchange coupling at the middle of the diamond is set 
to g{x = 0) = 0.5 in figures (a), (c), (e), and (g) and to 
g{x = 0) = 1.5 in (b), (d), (f), and (h). In (a), the YSR 
states are seen to cross zero energy when the effective ex¬ 
change coupling increases as x is changed from 0 to 1 towards 
the charge degeneracy points, (b) With a stronger coupling, 
the eye-shaped feature in the middle becomes smaller, and for 
strong enough coupling it disappears, (c) With finite phase 
difference, two new bound states appear, (d) A finite phase 
difference reduces the effective exchange coupling and restores 
the eye-shaped crossings. The detailed phase evolution is 
shown in (e) and (f) for a cut in gate voltage correspond¬ 
ing to X = 0.5. While panels (a)-(f) are for a symmetrically 
coupled junction (0 = 7r/4), (g) and (h) show the phase de¬ 
pendence for an asymmetrically coupled junction {0 = 7r/3). 
Note that the finite-bias degeneracy at 0 = tt is lifted by the 
asymmetry. 
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understood from the position of the YSR states at a; = 0 
(where W = 0) and (f) = 0: 


Ea- = -a 


1-g" 

1 + 52- 


(38) 


For (a), where ^ = 0.5atx = 0, the upper YSR state 
moves down and crosses zero as x approaches the charge 
degeneracy points at which g diverges. In contrast, in (b) 
where g = 1.5 at x = 0, the YSR state has already crossed 
zero at X = 0 due to stronger coupling. Interestingly, the 
zero energy crossings correspond to a change of parity 
of the ground state, cf. the discussion in Ref. [H Adding 
now a finite phase difference, panels (c) and (d) reveal the 
2-channel nature of the problem, with two bound states 
above and below the zero energy. 

In panels (e)-(h), we plot the dispersion of the bound- 
state with phase difference. In general, a finite phase 
difference is seen to to shift the value of the coupling at 
which the energy levels cross zero to higher values. In 
the plots, 0 = TT14 corresponds to symmetric coupling, 
th = see Eq. ([T^. For the particle-hole symmetric 
point X = 0 and symmetric coupling 6> = 7r/4, we see that 
the two excitations are degenerate at (j) = n. 

It is interesting to compare the above features with the 
experimental results by Chang et a/.^ They show con¬ 
ductance plots similar to Fig.[^(a)-(b) for three different 
ranges of back-gate voltage, corresponding to different 
devices with either strong, or weak coupling, as well as 
one device right at the transition where the YSR states 
touch at zero energy in the middle of the diamond at 
X = 0. Furthermore, Chang et al. show the phase depen¬ 
dence for different fixed gate voltages similar to Fig. 
(e) and (f) and with similar qualitative features: only 
weak gate dependence in the weakly coupled device, and 
a restoring of the zero-energy crossings at finite phase 
difference for the stronger coupled devices. 

We have seen that a finite phase difference results in 
two subgap states at positive energy. Interestingly, the 
same situation occurs for a magnetic impurity coupled 
to an s± superconductor, i.e., a superconductor with two 
bands, where the pairing potentials have different signs 
in the two bands.® Note that if the two superconduc¬ 
tors have different pairing potentials, but with same sign, 
there is only one bound state within the smallest of the 
two gaps. 


B. Conductance to the normal-metal tunnel probe 

The YSR subgap states derived above can be observed 
by means of tunneling spectroscopy from a normal metal 
lead, as illustrated in Fig. and experimentally carried 
out in Refs. and M We assume the Y-dot tunneling 
rate to be larger than any relaxation rate between the 
YSR state and the quasiparticle continuum in the super¬ 
conductors, thus ruling out the single-electron tunneling 
currents which were recently demonstrated to be impor¬ 
tant for analyzing scanning tunneling spectroscopy of Mn 



Fig. 4. (Color online) Differential conductance d//dF in 
units of 2e^//i between the normal electrode and the super¬ 
conducting leads, see Fig. The different lines correspond 
to constant phase cuts (0 = 0,7r/2,7r) in Fig. [^e)-(h), as in¬ 
dicated at the top of the figure. The left panels (a) and (c) 
are for weak coupling, (^(x = 0) = 0.5), while the right panels 
(b) and (d) are for stronger coupling, g(^x = 0) = 1.5. The top 
panels (a) and (b) are for symmetric coupling 0 = 7r/4, while 
the bottom panels (c) and (d) are for asymmetric coupling 
0 = 7T13. In all cases the gate voltage is set to the particle- 
hole symmetric point, x = 0, and the coupling to the normal 
lead was chosen to be gNN = 0.1. 


adatoms on a Pb(lll) surface.® In this case, the current 
is carried exclusively by elastic Andreev reflections and 
the differential conductance between the normal lead and 
the superconducting region can be obtained from the S- 
matrix in the normal lead a^®® (at zero temperature) 

^ = ^Tr[2-S,,Sl + S^Ae]^ (39) 

where See and She are the scattering amplitudes at en¬ 
ergy eV for an incident electron to be reflected as an elec¬ 
tron or a hole, respectively. The amplitudes See and She 
can be obtained from the retarded T-matrix in Eq. ( [27| ), 
which we write here as (setting 5 = 0 for simplicity) 

T^(u;) = V[1 - (40a) 

V = + Wr^ (40b) 


where r* and a* are Pauli matrices representing particle- 
hole and spin space, respectively, and J and W are 3x3 
coupling matrices in lead space, with elements Jaa' and 
Waa' for 0^5 = L, R, N. 

The unperturbed momentum-summed Green’s func¬ 
tion is a diagonal matrix in lead space, where for the 
superconducting leads a = L, R it is 

R (w + i7?) - A„(ricos(/;)-T2sin</>) 

+ -■ 

where 77 is a positive infinitesimal. For the normal lead, 
the Green’s function simplifies to = —iz^^Tr, 
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assuming the normal-lead density of states to be the 
same as for the two superconductors, since any difference 
can be absorbed into the tunneling matrix elements. The 
^'-matrix is now expressed through the T-matrix as 

S(cj) = 1 -27riz/i.T(cj), (41) 


and the amplitudes See and She are found as submatrices 
of S with a' = a = N and component (2,1) in electron- 
hole space for She and (1,1) for See- The calculation can 
be carried out analytically, but is quite lengthy. In Ap¬ 
pendix ^ we provide an analytical expression for dJ /dV 
when the coupling to the sup erconducting leads is sym¬ 
metric, 0 = 7r/4, see Eq. 

In Fig. we present the differential conductance when 
the subgap states are probed by the normal lead. The 
traces correspond to vertical cuts in Fig.|^ The coupling 
is chosen to be g^N = 0 . 1 : weak enough to resolve the 
YSR states as distinguishable conductance peaks, and 
large enough to actually see them. For the cases with two 
bound states, both states give rise to peaks in d//dF, 
but with different widths. Assuming (/) = 0,6> = 7r/4, 
X = 0, and gjsfN ^ Qll = dRR^ the width o f th e sub¬ 
gap conductance peak can be found from Eq. (Al) to be 
proportional to ^ArArA(l — F’^/A^), with Ecr from (36), 
implying very sharp peaks close to the gap-edges, and an 
overall scale set by the width ^fjVArA of a deep YSR state. 
When the voltage is resonant with subgap states the dif¬ 
ferential conductance is close to 2e^/h, except when the 
two subgap states are degenerate at zero enemy, in which 
case they add up to exactly 4e^/h (cf. Fig. ^(h)), or for 
a symmetric junction {0 = 7r/4) when they are degener¬ 
ate at finite energy for 0 = tt, where the conductance is 
exactly zero (cf. Fig. i(e-f)). 


C. Supercurrent 

At zero temperature, the supercurrent can be found 
as the derivative of the ground-state energy, Fgs, with 
respect to the phase difference between the two super¬ 
conducting leads 


/ 5 ( 0 ) = 2 


d(j) 


(42) 


If the ground-state energy has a single minimum either at 
0 = 0 or at 0 = TT, the junction is classified respectively as 
a 0-junction or a 7r-junction. In the former (0) case, the 
supercurrent is a continuous function of (j) with positive 
slope at 0 = 0 and negative slope at 0 = tt. In the 
latter (tt) case, the supercurrent has negative slope at 0 
and a positive slope at tt. In cases with minima both 
at 0 = 0 and at 0 = tt, a O'- or a 7r'-junction refers to 
the global minimum being respectively at 0 or tt, and 
the supercurrent is now a discontinuous function (with 
one discontinuity in the interval 0 G [0, tt] and one in the 
interval (j) G [7r,27r]) with positive slopes both at 0 = 0 
and at 0 = TT. 



Fig. 5. Phase diagram with couplings parameterized by 
Eq. (II), in terms of F = z^F(ti + ti^), and the dimension¬ 
less level position (gate voltage) x — 1 -\- 2edlU. The phase 
boundaries are independent of A, and their dependences on 
X and 0 are indicated, with the function f(0) defined as 
f{0) = [I/2{ sin^(2^) -h y^4 -h sin^(2^)}]^/^. In this diagram 
the coupling asymmetry parameter was chosen to 0 = tt/S, 
whereas for a symmetric junction with 0 = 7r/4, the O’-O- 
boundary will never be reached, consistent with earlier results 
on the Anderson model. ^ * ^° * ^^ * 



Fig. 6. (Color online) Supercurrent vs. phase difference at the 
particle-hole symmetric point, x = 0, for exemplary values of 
g and 0 (see inset) corresponding to the four different kinds 
of Josephson junctions. Similar results were first obtained in 
Refs. [23 and [29| and experimentally demonstrated in Ref. (52] 


Within the polarized-spin approximation, we find the 
continuum to be independent of the phase difference in 
the limit of infinite bandwidth. This is is akin to a so- 
called 5/z6>rt junctiorP^and implies that the supercurrent 
can be obtained from the subgap excitation spectrum 
alone, which according to the standard Bogoliubov-de 

























Gennes formalism is given by 


^GS — (l^+l ^ l^-l) + const. (43) 

The resulting phase diagram as a function of the dimen¬ 
sionless gate voltage x, and the lead-dot tunneling rate, 
r = is depicted in Fig. and examples of 

supercurrent in the different phases are shown in Fig. 
Only for small coupling constants do we have a sinusoidal 
current-phase relations: 

Is ~ (A/ 2 ) sin^(26>) — g^) sin 0 . (44) 

In some cases the supercurrent is discontinuous, which 
is directly related to the subgap states crossing zero en¬ 
ergy and changing their spin (as is evident from Eq. ([4^). 
Similar behavior of the supercurrent has been demon¬ 
strated experimentally in Ref. [52] where a thermally 
smoothened discontinuity was observed at a phase dif¬ 
ference (j)c whose dependence on gate gate voltage was 
shown to comply well with the formula 

(j)c = 2arccos[A/7 — (x/h)^]. (45) 


This functional dependence on the dimensionless gate 
voltage, X, was derived in the atomic limit (A ^ oo) in 
Ref. |30l However, in the experiment A ^ T^, and 7 and 
h were therefore used as free fitting parameters. More¬ 
over, the data was also shown to match the result of a 
quantum Monte Carlo calculation for the corresponding 
Anderson model. We note here that the same functional 


behavior follows directly from Eqs. (34), (42), and (43), 
with coefficients 7 and h defined as 


7 = 


1 


1 


tan^(26>)’ 


h = g sin( 2 ^). 


(46) 


which allows extracting the coupling asymmetry and the 
dimensionless exchange coupling. 


the excited subgap singlet perturbatively matches the 
second order expansion (in gaa') of the excitation en¬ 
ergy obtained from the polarized-spin approximation for 
H = 0 and S = when replacing ^ ‘^g^cx' in (34). 


In Appendix |B 2 [ we show how the subgap excitation 
energies at weak coupling, g are shifted in energy 

by an external magnetic field. Assuming that the field 
is screened in the superconducting leads, g^e ~ 0 , the 
transition energy from the ground-state doublet to the 
excited-state singlet is found to approach the polarized- 
spin result in the limit of high magnetic field, B ^ ^^|A| 
(seeEig. [^. 

Appendix |B 3| goes beyond the low-energy single¬ 
quasiparticle Hamiltonian and provides a calculation of 
the leading order correction to the ground-state 

energy. With reference to the three-quasiparticle Yosida 
wavefunction ansatz analyzed in Ref. |38l we argue why 
this shift of ground-state energy does not modify the ex¬ 
citation energy found in Appendix |B 1 | 

Einally, this second order shift in ground-state energy 
is used in Appendix |B 4| to calculate the supercurrent in 
the presence of a finite magentic field via formula (42). In 
the perturbative regime, the functional dependence of the 
supercurrent on both (j) and 0 is similar with, and without 
the polarized-spin approximation, but again, only in the 
large-held limit do they match exactly. 

Summarizing theses hndings, we hnd qualitative agree¬ 
ment between the results found within the polarized-spin 
approximation and the perturbative results for the full 
quantum mechanical treatment. A non-perturbative cal¬ 
culation of the YSR spectrum and the supercurrent be¬ 
yond the polarized-spin approximation could be found by 
means of a numerical renormalization group calculation 
of the T-matrix, whereas the non-linear Andreev conduc¬ 
tance would require an average of the current operator. 


V. CONCLUSIONS 


IV. QUANTUM MECHANICAL TREATMENT 
OF THE SPIN DEGREE OF FREEDOM 

Here we discuss the results for the fully quantum me¬ 
chanical description, where unlike for the polarized-spin 
approximation an exact solution is not possible. There¬ 
fore we resort to perturbation in the cotunneling cou¬ 
plings g and re, which is valid when the Kondo tem¬ 
perature is much smaller than the superconducting gap 
Tk ^ A and the ground state is a doublet. Details of 
the calculations can be found in Appendixand here we 
simply summarize the main findings and contrast them 
to the results for the polarized spin found above. 

In Appendix |B 1 [ we diagonalize an effective low- 
energy single-quasiparticle Hamiltonian. As in the spin- 
polarized case, we find two subgap excitation energies 
when the two superconductors have a phase difference. 
The transition energy from the ground-state doublet to 


In summary, we have determined the phase-dispersion 
of the Yu-Shiba-Rusinov states induced by a spinful 
Coulomb blockaded quantum dot coupled to two phase- 
biased superconducting leads. At finite phase difference, 
two channels are involved in the screening of the dot 
spin. Consequently, the phase-biased system exhibits 
two, instead of one, YSR states, one of which merges 
with the continuum at the gap edge for zero phase dif¬ 
ference. We have shown how the corresponding subgap 
excitation spectrum is modified by coupling asymmetry 
and potential scattering, and established that the phase 
difference generally shifts the parity transition, and the 
accompanying sign change in the supercurrent, to larger 
values of the exchange coupling. 

We have solved the problem exactly in the spin- 
polarized approximation, and perturbatively in the fluc¬ 
tuating quantum case (cf. Appendix]^, and obtained 
a closed analytical expression for the subgap excitation 
energy [Eq. ([M|)], which extends Rusinov’s resulfP to a 
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situation with two superconductors at a finite phase dif¬ 
ference and arbitrary tunnel couplings. 

The YSR bound states can for example be observed by 
a normal metal tunnel probe connected to the dot, which 
would also give information about their spectral weight. 
We have therefore calculated the differential conductance 
in such a setup. These results should provide a valuable 
basis for a more detailed analysis of future experiments 
like Refs. [ 12 ] and [T3l 
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Appendix A: Explicit formula for the 
symmetric-coupling conductance 

In the case of symmetric coupling of the quantum dot 
to the two superconductors {0 = 7 r/ 4 ), we can derive a 
simple closed-form expression for the differential conduc¬ 
tance d//dU. For subgap conductance (|U| < A) we 
obtain 


dl 4t+t_/^cos^ ^ X 2e^/h 

dV [1 + - r+r_)(l - p sin^ |) + (r+ -h r_)/u]2 + PP- + “ P sin^ |) - (t+ + tpfv]‘^ 

and for continuum conductance (|U| > A) we get 


+(u ^ -u), 
(Ala) 


dU 


2f\v\[{r^t_ + rH+){l + p sin^ |) + (4 + ^_){1 + + p sin^ f)} + 4:t+t_f\v\] x 2e^//i 

[1 + {t+t- - r+r_)(l + /2 sin^ f) + {t+ + t-)f\v\]‘^ + [{r+t- + r_i+)(l + p sin^ |) + (r+ + r_)/|i;|]2 

(Alb) 


with V = V/A and f{V) = A/\/|A^ — and where 

transmission, and reflection amplitudes for incident elec¬ 
trons in the superconductors have been identified as 

. _ ‘^{dNS Y Wns)‘^ / * 

^ l^{gNN±WNNy^ ^ ^ 

r± = g ±w — {gNN ± WNN)t±. (A3) 


Matsuura and Nagaoka,!^ using Yosida’s wave-function 
ansatz.l^ Since the calculation using the ansatz is rather 
technical, it is more instructive to calculate the subgap 
spectrum, using an effective single-quasiparticle model 
for the case with no potential scattering at x = 0 . 

First, the BCS leads are diagonalized by the Bogoli- 
ubov transformation 


where gNS = gNL = gNR and wns = ^nl = wnr- 
Similar expressions to Eq. were derived for polar¬ 

ized spin coupled to a single superconductoit^ and un¬ 
conventional superconductor junctions containing sub¬ 
gap states.l^ For calculations of d//dU in Fig. 3 when 
there is finite coupling asymmetry, the couplings involv¬ 
ing the normal lead are parametrized as 


gNL = \/gNN g cosO, 
WNa = gNa'^tg, 


gNR = V^ArAr^sin6>, 
"^NN = gNN'^/g- 


Appendix B: Subgap states beyond the 
polarized-spin approximation 

In this appendix, we determine the subgap states in¬ 
duced by the dot spin without doing the polarized-spin 
approximation applied in Section |in| The presence of the 
spin-flip terms prohibits a solution along the lines in the 
main part of the paper and we therefore revert to a per¬ 
turbative treatment, valid for small dimensionless cou¬ 
plings. This problem was considered already by Soda, 


Cafc<T = (Bl) 

where 



(B2) 

In terms of the Bogoliubov quasiparticle operators 
the lead Hamiltonian <§ reads 

■^LR ^ V '^o^fecrTg/ccr'Ta/ccr; 
aka 

with eigenenergies 

+ = + (B4) 


1. Effective single-quasiparticle model 

Next, we express the exchange Hamiltonian Hj in 
terms of the Bogoliubov operators. In accordance with 
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the leading order term in Yosida’s ansatz, we neglect all 
pairing-like terms, and in the exchange 

Hamiltonian and obtain the following low-energy effec¬ 
tive model: 


Hj 


(B5) 


a k'ak 

^ a'k'akx) 

+ + S~llk'flaki 

= '^i^lM'ipk, 

k'k 


where we have set ^ ^ ^ since their energy depen¬ 

dence only matters for higher order corrections to the 
subgap excitation energies. The last line in (B5) is ex¬ 


pressed in terms of the conduction electron 4-spinor 
~ ('yLkX^"^Rkt^"^LkV"fRki)^ 

and the matrix 

^ S^\^f Jr.r. 

7 - -^v 

= Ms (8) At/, 


M = 


Jll 

7* l+e“ 
y^LR 2 


LR- 


Jrr 


(B7) 


in which Ms operates in spin, and Mi in lead space. 
The Hamiltonian (B5) is written in the excitation ba¬ 


sis with 7 ^^^ annihilating the BCS vacuum, 7Q,/eo-|0) = O 5 
and it can be diagonalized exactly. After diagonalizing 
the lead space matrix Mi we obtain two decoupled chan¬ 
nels S = 1,2 


where 


7l/ca ~ ^^L/ccr + 
"f2ka = ^7i?/ccr ~ 




7 ^ f -i Jll — Jrr 

‘=vn ‘—— 


(B 8 a) 


(B 8 b) 


are the singlet 

7^*;) = ^ (|t<5fe,4-) - |4'<5fe,t)) 5 V =- 2 ) (BlOa) 
and triplet states 


l7fc) = (lt<5fe,4-) + |4'<5fe,t)) , 

\TXk) = nskA), \Ts-k) = \Uk,i), At = ^. 


(BlOb) 


Here \s/t denotes the corresponding eigenvalues and 

=7lfea|0)l7, with ^^Is) = s/2|s). (Bll) 


We note that the states dBl^ and ( |B11| ) span only the 
Hilbert sub-space of single particle excitations with re¬ 
spect to the ground-state doublet \ Ds) = |0)|s). Ex¬ 
pressed in the singlet/triplet basis ( B1Q[ ), the Hamilto¬ 
nian now takes the following simple form 

ifLR + Hj = J2Ek {\Ssk){Ssk\ + 

6k 

-lT.L\Ssk'){Ssk\ + lYl 7|T7)(Ti,| 


Sk'k 


jSk'k 


(B12) 

To fi nd a singlet subgap state from the above Hamil¬ 
tonian (B12) we form the linear superposition 


is<> = E ^6k\Ssk), 


(B13) 


and solve the stationary Schrodinger equation 

{Hi^r^Hj-E)\Ss)=0. (B14) 


Projecting Eq. (B14) to {Ssq\ we obtain the equation 

3 J(5 Ask 


^5q — 


^ 2 Eq-E' 


which is integrated over q to yield 
1 = 


(B15) 


(B16) 


The corresponding eigenvalues for two channels are 


Ji/2 = X {Jll + Jrr ± Jd ), 

I - 6 

Jd = y {Jll - Jrr)^ + ^Jlr\^ cos^ -. 


The eigenstates of the spin matrix 

|t(5/c,t) |t(5/c,i) |i(5/c,t) 

( \ 0 0 

0 1 

0 1 -1 

\ 0 0 0 


11(5/e, I) 


0 

0 

0 

1 

2 


\ 

, (B9) 

/ 


where gs = evfJsI2 and the necessary integral, Ie^ for 
subgap states with |E| < A and large bandwidth D ^ A 
is given by 


Ie 


1 

q 





^ ‘^E{\ + Aarcsinf) 

VA2 -7;2 


(B17) 


We note that the subgap tripl et solut ions can exist for 
ferromagnetic coupling, gs < However, the An¬ 

derson model always gives rise to antiferromagnetic ex¬ 
change, and therefore there will be no triplet subgap 
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state. Parameterizing the leading-order energy of singlet 
solutions as 


Eos = A (1 - (B18) 

we finally obtain the perturbative solution valid to lowest 
(second) order in gs 


We start by examining Eq. ( |B25 ) in the low magnetic 
field limit where \B\ Is'lAj. Parameterizing the energy 
as before, 


— A(1 — 77I), 

we first expand Eq. ( B25[ ) to lowest order ir 


(B26) 


\r]os\ = ^V2gs. (B19) 


This result matches the second order expansion (in gaa') 
of the excitation energy obtained from the polarized-spin 
approximation for re = 0 , 5 = 0 , and S' = |, when the 
replacement g^a' ^9aa' is made in (34). This could be 
anticipated already by comparing the expectation value 
of Hj for the state |t( 5 /c:i) corresponding to a polarized 
spin and for the state \Ssk)'- 


{S6k\Hj\Ssk) 


(B20) 


»7± = y»?|±f, (B27) 

to obtain the equation 

-\/2g5(7?+ + r?_)-6fit| =0. (B28) 

Expanding now to lowest order in this yields 

rfs - - ^gWs - ( 7^1 = 0, (B29) 


If the potential scattering term, Hw^ is included, the 
perturbative result and the polarized-spin approximation 
still match (see Ref. l38|) . 


2. Magnetic field dependence 


In the singlet/triplet basis, the Zeeman term takes the 
following form 

Hb = -BY,{\Ssk){T^k\ + \T!k){Ssk\) 

6k 

+ BY.iKk)iTs\\ - |r,-,)(r,-,|), (B21) 

6k 

where the last term represents the Zeeman splitting of 
the ground-state doublet |T)s) = |0)|s). In terms of the 
two potentially different ^-factors, we have introduced 
the difference and average 5-fields as 

B=^{gd-gce), B=^{gd + gce)- (B22) 


Eor gd ^ we have 5 7 ^ 0 and the singlet and the 
triplet are mixed to form a new eigenstate 

IV-^) = + bsk\Tl)) . (B23) 

k 

Projecting the stationary Schrodinger equation 

{Hbr + Hb+Hj- E)\i;s) = 0 (B24) 

to the sub-space spanned by \Ssq) and \T^q) and integrat¬ 
ing over q yields the secular equation 


^ 2^ + ^E-b\ 

^[^e+b ~ ^e-b\ 


^2 \.^E+B ^E-b] 

I -h ^[^e+B + ^E-b] 


= 0 . 
(B25) 


which has the following leading order perturbative solu¬ 
tion 


+ (B30) 

This shows that the energy of the subgap state decreases 
quadratically with B for small magnetic fields. 

In the h igh field limit where |5| ^ |^|A|, we find from 
Eq. ( B25[ ) to lowest order in gs that 


Es = A{1B - B, \gs,c]\^ V2gs, 


(B31) 


which correspon ds to neglecting off-diagonal terms in 
the spin matrix (B9). Eor intermediate magnetic field 


strengths, the perturbative (in g s) solution is obtained 
by numerically solving Eq. (B28). 


Einally, we examine the subgap excitation spectrum in 
the case where the magnetic field in the superconducting 
leads is screened, i.e., for gee ~ 0 and B = B = gdBj2. 
Note that we assume the magnetic field to be much 
weaker than the critical field and hence neglect its in¬ 
fluence on the gap. Eor positive magnetic field, 5 > 0, 
the ground state is the lower energy component of the 
doublet with 5^ = — 5 and the resulting subgap exci¬ 
tation energy, — 5^, for a particular channel 

S is depicted in Eig. Eor high magnetic fields, this 

excitation approaches the energy of the polarized spin 
approximation, Eq. ( [M| ). In this case the eigenstate is 

|t<54) = Efc ^Llt<5fe4), where is determined by pro- 
jecting the Schodinger equation to |t( 5 /c, i) with neglected 
spin -flip terms S~^ and S~ in the exchange Hamiltonian 
(B5). Changing the sign of the magnetic held simply 


reverses all spins in the previous discussion. 


3. Beyond the single-quasiparticle approximation 

So far we have examined the excitation spectrum when 
only a single quasiparticle is included and the system is 
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Fig. 7. (Color online) Subgap excitation spectrum, = 
Es — with respect to the ground-state doublet \Di) for 
particular channel (5 = 1,2. It is assumed that the magnetic 
field in the superconductor is screened pee ~ 0. For antifer¬ 
romagnetic coupling, gs = 0.05, the excited state is singlet 
like, which for high magnetic fields becomes polarized-spin 
like |t5 4)- 


effectively described by the Hamiltonian (B5). We now 
return to the original Hamiltonian and investigate the ef¬ 
fects of the terms 7^/^7^^ which were neglected to arrive 
at (B5). With these terms included, second order per¬ 
turbation theory yields the following shift of the ground- 
state energy from which the supercurrent can be deduced 


7^^ = E 


\H'\D, 


7 °^ - El 

1 Jcc>ceJ‘^ak‘^a'k'^ 


■E 

ak 

a'k' 


E. 


ak 


' ^a'k' 


^ 1 Ja'aJaa'l'^a'k'^^ak^ ^^ ^ak'^a'k' ^ \ 

2 + Eq,/]^/ s(^g g^o^B 

^ ^^a'a^^aa'l^^a'k'"^ak^ E Uak'^a'k'\ \ 

Eak E E(yfj^f J 

(B32) 

with |/) denoting all possible intermediate states and 

H' = Hj^ Hw. 

This second order shift in ground-state energy would 
appear to influence the observable excitation energies, 
but when analyzed in terms of a three quasiparticle 
Yosida wavefunction ansatz (see Ref. [38]), the eigenen- 


ergies of the subgap states are found to b e shi fted in 
exactly the same way, and overall the shift (B32) drops 
out and the previously obtained second order result for 
the excitation energies, Eq. (B19), remains valid. In this 


manner, the Yosida wave-function ansatz generates a well 
defined perturbative expansion for the energy differences 
in the dimensionless couplings, g and w. 

4. Supercurrent 


Going beyond the polarized-spin approximation, we 
calculate the supercurrent perturbatively from Eqs. (42) 


and (B32). Eor H > 0 the ground state is \D^) and one 
finds 


/^ = 2A sin^(26>) — p^) F(0) — 2p^F(H) sin0. 


to leading order in w and p and with 

1 UkVkUk'Vk' 

^ ~ A_2,.2 


Ek E Ef + B 


(B33) 

(B34) 


whereby E(0) ~ 1/4 for D ^ A. Eor zero magnetic field, 
the supercurrent simplifies to 

^ sin^(26>) (w‘^ — 3p^) sin0, (B35) 

which always corresponds to a 7r-junction since g > w. 
Eor small magnetic fields with \B\ A, we have 


F{B)-F{0)^-^. (B36) 

Erom Eq. ( |B33[ ), it then follows that a positive B de¬ 
creases, and a negative B increases the magnitude of the 
supercurrent, \Is\- Notice that B can become negative 
for positive B when pee > Qd- 

Eor large magnetic fields with \B\ ^ A, one finds in¬ 
stead 


E{B) 


A 

B 




(B37) 


and as E{B) vanishes with increasing field, the spin is 
polarized, and the supercurrent (44), obtained within the 
polarized-spin approximation, is recovered. 
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